########################################################################
# PLS 206 — Week 2 live coding in class
# Grey Monroe, Fall 2026
########################################################################

########################################################################
# 1. Read the data and plot it first
########################################################################

# A winter cover crop is sown in fall and ended in spring; the more
# biomass it makes by spring, the more it protects the soil and feeds
# the next crop. cover_crop.csv: spring dry biomass (kg/ha) of the same
# cover crop on 50 farms (simulated), with each farm's
#   location:  latitude, longitude, elevation (m)
#   winter weather, Oct-Mar: gdd (growing degree days, base 4 C),
#     tmean_winter (C), tmin_coldest (min temp of coldest month, C),
#     frost_days, rain_total (mm), rain_wettest, rain_driest (mm in the
#     wettest / driest month), dry_spell (longest run of dry days)
#   soil:      sand, silt, clay (%), soil_ph, organic_matter (%)

# install.packages(c("ggplot2", "car"))
farms <- read.csv("https://greymonroe.github.io/PLS_206/data/cover_crop.csv")

str(farms)
mean(farms$biomass)
sd(farms$biomass)
hist(farms$biomass)
plot(farms$gdd, farms$biomass)

########################################################################
# 2. Fit a line and read three numbers
########################################################################

# In summary(), the slope is the Estimate on the gdd row (kg/ha of
# biomass per growing degree day), its p-value is Pr(>|t|) on the same
# row, and R² is Multiple R-squared near the bottom. confint() gives
# the 95% confidence interval for the slope.

model <- lm(biomass ~ gdd, data = farms)
summary(model)
confint(model)

# The same numbers pulled out one at a time:

coef(model)               # intercept, then slope
confint(model)["gdd", ]   # CI for the slope
summary(model)$r.squared
sigma(model)              # residual standard error: the typical miss, in kg/ha

# Does winter rainfall predict biomass?

model_rain <- lm(biomass ~ rain_total, data = farms)
summary(model_rain)
confint(model_rain)

########################################################################
# 3. A correlation test is the same test
########################################################################

# With one predictor, cor.test() gives the same t and p as the slope,
# and r² is R².

cor.test(farms$gdd, farms$biomass)
cor(farms$gdd, farms$biomass)^2

########################################################################
# 4. Draw the line
########################################################################

plot(farms$gdd, farms$biomass, xlab = "Growing degree days", ylab = "Biomass (kg/ha)")
abline(model, col = "red", lwd = 2)

library(ggplot2)
ggplot(farms, aes(gdd, biomass)) +
  geom_point() +
  geom_smooth(method = "lm") +
  theme_classic()

########################################################################
# 5. Check the fit
########################################################################

# Residuals vs fitted: no curve, even spread. Q-Q: points on the line.
# Residuals vs leverage: nothing outside the dashed Cook's distance
# lines.

par(mfrow = c(2, 2))
plot(model)
par(mfrow = c(1, 1))
shapiro.test(resid(model))

########################################################################
# 6. Predict
########################################################################

# newdata is a data frame whose column has the same name as the
# predictor.

predict(model, newdata = data.frame(gdd = 1000))
range(farms$gdd)                                  # where the data are
predict(model, newdata = data.frame(gdd = 0))     # outside the data: don't trust it

# Two intervals around the same prediction (fit):
# confidence: where the LINE is, i.e. the mean biomass of all farms with
#   1000 GDD. Only uncertainty in the slope and intercept, so it is narrow,
#   and it shrinks with more data.
# prediction: where ONE NEW FARM with 1000 GDD will land. The line's
#   uncertainty plus that farm's own scatter around the line (sigma), so
#   it is much wider, and more data never shrinks it below about +/- 2 sigma.

new <- data.frame(gdd = c(800, 1000, 1200))
predict(model, newdata = new, interval = "confidence")
predict(model, newdata = new, interval = "prediction")

# A warmer winter: every farm gets 100 more growing degree days.
# The mean prediction moves by exactly 100 x the slope.

warmer <- predict(model, newdata = data.frame(gdd = farms$gdd + 100))
mean(warmer)
mean(warmer) - mean(farms$biomass)
100 * coef(model)["gdd"]

########################################################################
# 7. Screen many predictors at once
########################################################################

# One correlation per column, without copying cor() over and over.

predictors <- setdiff(names(farms), c("farm", "biomass"))
cors <- sapply(farms[, predictors], function(v) cor(v, farms$biomass))
round(sort(cors), 2)
names(which.max(abs(cors)))   # the strongest, ignoring sign

# the same with a for() loop
cors <- c()
for (p in predictors) cors[p] <- cor(farms[[p]], farms$biomass)

########################################################################
# 8. Two predictors
########################################################################

# Farms farther north get fewer growing degree days. Is it latitude
# itself, or the warmth that goes with it? In a model with both, each
# slope is the effect of that predictor holding the other fixed. vif()
# checks whether the two overlap too much to separate (rule of thumb:
# above 10 is a problem).

cor(farms$gdd, farms$latitude)
summary(lm(biomass ~ latitude, data = farms))
model2 <- lm(biomass ~ latitude + gdd, data = farms)
summary(model2)

library(car)
vif(model2)
anova(model, model2)   # does latitude improve on gdd alone?

# Higher farms also have colder winters. Elevation alone, then with gdd:

summary(lm(biomass ~ elevation, data = farms))
summary(lm(biomass ~ elevation + gdd, data = farms))

########################################################################
# 9. Interactions
########################################################################

# Does warmth help more on wetter farms? gdd * rain_total fits both
# slopes plus their product; the gdd:rain_total row tests whether the
# gdd slope changes with rainfall.

model_int <- lm(biomass ~ gdd * rain_total, data = farms)
round(coef(summary(model_int)), 4)
# no interaction here: keep the simpler model

########################################################################
# 10. Every predictor at once
########################################################################

# ~ . means every other column; - farm leaves out the ID column. One
# coefficient comes back NA: look at the soil columns and think about
# what sand, silt, and clay have to add up to.

full <- lm(biomass ~ . - farm, data = farms)
summary(full)
summary(farms$sand + farms$silt + farms$clay)

# Drop the redundant column and compare with the one-predictor model.
# R² can only go up as you add predictors, so compare with adjusted R²,
# and check vif().

full2 <- lm(biomass ~ . - farm - clay, data = farms)
summary(full2)
summary(model)$adj.r.squared   # gdd alone
summary(full2)$adj.r.squared   # everything
round(vif(full2), 1)

########################################################################
# 11. Which predictor matters more?
########################################################################

# Slopes in different units can't be compared. scale() puts every
# variable in standard deviations, so the slopes become comparable.

model_std <- lm(scale(biomass) ~ scale(latitude) + scale(gdd), data = farms)
round(coef(model_std), 2)

########################################################################
# 12. Report it
########################################################################

# An effect size in units with its CI, the test, R², and n: the numbers
# for one sentence in a results section.

coef(model)["gdd"]
confint(model)["gdd", ]
coef(summary(model))["gdd", ]   # t and p
summary(model)$r.squared
nobs(model)